#include "gpuSpokeDart.h"
#include "cuda_defs.h"
#include "curand_kernel.h"
#include "randUtilGpu.h"
#include "randUtil.h"
#include "device_launch_parameters.h"

const unsigned int MAX_GRID_CELL_ELEMENTS = 8;

namespace GPUSpokeDartCuda
{
__device__ __constant__ DATA_TYPE c_radius;
__device__ __constant__ DATA_TYPE c_sqRadius;
__device__ __constant__ unsigned int c_m;
__device__ __constant__ DATA_TYPE c_invNumDimensions;

__device__ unsigned int d_checkedSamples[NUM_BLOCKS];
__device__ unsigned int d_front[MAX_FRONT];
__device__ DATA_TYPE2 d_intersectionArray[NUM_BLOCKS*NUM_THREADS*MAX_INTERVAL];
__device__ unsigned int d_numOverlaps[NUM_BLOCKS*NUM_THREADS];
__device__ int d_lock = 0;
__device__ unsigned char d_hasSolution[NUM_BLOCKS*NUM_THREADS];

// for neighbor query (2D only)
__device__ __constant__ DATA_TYPE d_cellSize;
__device__ __constant__ int d_gridMax; // ceil(1.0 / cellSize)
unsigned int* d_numGridCellElements;
unsigned int* d_gridcells;
__device__ void getGridMinMax(DATA_TYPE point, int& minIndex, int& maxIndex)
{
	minIndex = (int)(point / d_cellSize) - 2;
	maxIndex = minIndex + 5;
	minIndex = max(0, minIndex);
	maxIndex = min(d_gridMax, maxIndex);
}
__device__ void addToGrid(VectorGPU& sample, unsigned int index, unsigned int* cells, unsigned int* numCellElements)
{
	int x = (int)(sample[0] / d_cellSize);
	int y = (int)(sample[1] / d_cellSize);
	int pos = (y * d_gridMax + x);
	unsigned int nextIndex = numCellElements[pos];
	cells[pos * MAX_GRID_CELL_ELEMENTS + nextIndex] = index;
	numCellElements[pos] = nextIndex + 1;
}
////////////////////////////////////////

__global__ void initRandomStates(curandState* randStates, unsigned int* pNextFront, GPUDart::GPUMDVector samples, VectorGPU* pFirstSample)
{
	unsigned int id = threadIdx.x + blockIdx.x * blockDim.x;

	curandState* localState = randStates + id;
	curand_init(RANDOM_SEED, id, 0, localState);

	if (id == 0)
	{
		d_front[0] = 0;
		*pNextFront = 1;

		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			samples(0, d) = (*pFirstSample)[d];
		//printf("First sample : %f %f\n", samples(0, 0), samples(0, 1));
	}
}

__global__ void setupGPUDarts(unsigned int* pNumFronts, GPUDart::GPUMDVector samples, curandState* randStates, GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions,
	GPUDart::GPUInterval* intervals)
{
	unsigned int numFronts = *pNumFronts;
	if (numFronts == 0)
		return;
	unsigned int id = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int trialId = id / numFronts;
	unsigned int frontId = id - trialId * numFronts;
	
	curandState* localState = randStates + id;

	DATA_TYPE refPoint[DIMENSIONS];
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		refPoint[d] = samples(d_front[frontId], d);

	if (trialId == 0)
	{
		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			refPoints(frontId, d) = refPoint[d];
	}

	VectorGPU direction;
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
	{
		direction[d] = RandUtilGpu::randGaussian(localState);
	}
	direction.normalize();
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		directions(id, d) = direction[d];

	GPUSpokeDart::GPUInterval& interval = intervals[id];
	DATA_TYPE intervalMin = -2 * c_radius;
	DATA_TYPE intervalMax = 2 * c_radius;
	for (unsigned int i = 0; i < DIMENSIONS; ++i)
	{
		if (direction[i] > 0)
		{
			DATA_TYPE t = ((DATA_TYPE)1 - refPoint[i]) / direction[i];
			intervalMax = min(intervalMax, t);

			t = (0 - refPoint[i]) / direction[i];
			intervalMin = max(intervalMin, t);
		}
		else if (direction[i] < 0)
		{
			DATA_TYPE t = (0 - refPoint[i]) / direction[i];
			intervalMax = min(intervalMax, t);

			t = ((DATA_TYPE)1 - refPoint[i]) / direction[i];
			intervalMin = max(intervalMin, t);
		}
	}
	interval.prev[0] = -1;
	interval.head = 0;
	interval.tail = 0;
	interval.numIntervals = 0;
	interval.numUsed = 1;
	if (intervalMin < -c_radius)
		interval.insert(interval.tail, intervalMin, -c_radius);
	if (intervalMax > c_radius)
		interval.insert(interval.tail, c_radius, intervalMax);
}

__global__ void updateIntervalsKernel(unsigned int* pNumFronts, GPUDart::GPUMDVector samples, GPUSpokeDart::GPUInterval* intervals,
	unsigned int *d_numSamples, GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions,
	unsigned int* d_gridcells, unsigned int* d_numGridCellElements)
{
	unsigned int numFronts = *pNumFronts;
	if (numFronts == 0)
		return;
	unsigned int id = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int frontId = id - (id / numFronts) * numFronts;

	unsigned int numSamples = *d_numSamples;

	if (threadIdx.x == 0)
		d_checkedSamples[blockIdx.x] = numSamples;

	DATA_TYPE refPoint[DIMENSIONS];
	DATA_TYPE direction[DIMENSIONS];

	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		refPoint[d] = refPoints(frontId, d);
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		direction[d] = directions(id, d);
	
	DATA_TYPE2* __restrict__ intersectionArray = d_intersectionArray + id;
	unsigned int numOverlaps = 0;

#ifdef USE_GRID_BASED_NEIGHBOR_SEARCH
	// 2d only

	int xMin, xMax, yMin, yMax;
	getGridMinMax(refPoint[0], xMin, xMax);
	getGridMinMax(refPoint[1], yMin, yMax);
	//printf("[%d] %d check with [%d-%d][%d-%d]\n", id, frontId, xMin, xMax, yMin, yMax);
	DATA_TYPE otherSample[DIMENSIONS];
	for (int y = yMin; y < yMax; ++y)
	{
		for (int x = xMin; x < xMax; ++x)
		{
			int pos = y * d_gridMax + x;
			unsigned int numElements = d_numGridCellElements[pos];
			for (unsigned int i = 0; i < numElements; ++i)
			{
				unsigned int otherSampleIndex = d_gridcells[pos * MAX_GRID_CELL_ELEMENTS + i];
				{
					//printf("[%d] %d check with %d\n", id, frontId, otherSampleIndex);
					for (unsigned int d = 0; d < DIMENSIONS; ++d)
						otherSample[d] = samples(otherSampleIndex, d);
					// ray-sphere intersection
					// Compute A, B and C coefficients
					DATA_TYPE oc[DIMENSIONS];
					for (unsigned int d = 0; d < DIMENSIONS; ++d)
						oc[d] = refPoint[d] - otherSample[d];

					DATA_TYPE a = 0;
					for (unsigned int d = 0; d < DIMENSIONS; ++d)
						a += direction[d] * direction[d];

					DATA_TYPE b = 0;
					for (unsigned int d = 0; d < DIMENSIONS; ++d)
						b += oc[d] * direction[d];
					b *= 2.0;

					DATA_TYPE c = -c_sqRadius;
					for (unsigned int d = 0; d < DIMENSIONS; ++d)
						c += oc[d] * oc[d];

					//Find discriminant
					DATA_TYPE disc = b * b - 4 * a * c;

					// if discriminant is negative there are no real roots, so return 
					// false as ray misses sphere
					if (disc >= 0)
					{
						intersectionArray[numOverlaps * NUM_THREADS * NUM_BLOCKS].x = b;
						intersectionArray[numOverlaps * NUM_THREADS * NUM_BLOCKS].y = c;
						++numOverlaps;
					}
				}
			}
		}
	}
#else

	__shared__ DATA_TYPE sharedSamples[DIMENSIONS][NUM_THREADS];
	DATA_TYPE localCache[DIMENSIONS];
	
	// load samples to register
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		localCache[d] = samples(threadIdx.x, d);
	
	for (unsigned int j = 0; j < numSamples; )
	{
		// copy register to shared memory
		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			sharedSamples[d][threadIdx.x] = localCache[d];

		__syncthreads();

		j += NUM_THREADS;

		// load samples to register
		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			localCache[d] = samples(j + threadIdx.x, d);

		int numIterationSamples = min(NUM_THREADS, numSamples - (j - NUM_THREADS));
		
		for (unsigned int i = 0; i < numIterationSamples; ++i)
		{
			// ray-sphere intersection
			// Compute A, B and C coefficients
			DATA_TYPE oc[DIMENSIONS];
			for (unsigned int d = 0; d < DIMENSIONS; ++d)
				oc[d] = refPoint[d] - sharedSamples[d][i];

			DATA_TYPE a = 0;
			for (unsigned int d = 0; d < DIMENSIONS; ++d)
				a += direction[d] * direction[d];

			DATA_TYPE b = 0;
			for (unsigned int d = 0; d < DIMENSIONS; ++d)
				b += oc[d] * direction[d];
			b *= 2;

			DATA_TYPE c = -c_sqRadius;
			for (unsigned int d = 0; d < DIMENSIONS; ++d)
				c += oc[d] * oc[d];

			//Find discriminant
			DATA_TYPE disc = b * b - 4 * a * c;

			// if discriminant is negative there are no real roots, so return 
			// false as ray misses sphere
			if (disc >= 0)
			{
				intersectionArray[numOverlaps * NUM_THREADS * NUM_BLOCKS].x = b;
				intersectionArray[numOverlaps * NUM_THREADS * NUM_BLOCKS].y = c;
				++numOverlaps;
			}
		}
		__syncthreads();
	}
#endif
	d_numOverlaps[id] = numOverlaps;
}

__global__ void updateIntervals2Kernel(unsigned int* pNumFronts, GPUDart::GPUMDVector samples, GPUSpokeDart::GPUInterval* intervals,
	unsigned int *d_numSamples, GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions,
	unsigned int* d_gridcells, unsigned int* d_numGridCellElements)
{
	if (*pNumFronts == 0)
		return;
	unsigned int id = threadIdx.x + blockIdx.x * blockDim.x;

	GPUDart::GPUInterval& interval = intervals[id];

	DATA_TYPE2* __restrict__ intersectionArray = d_intersectionArray + id;
	unsigned int numOverlaps = d_numOverlaps[id];

	DATA_TYPE direction[DIMENSIONS];
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		direction[d] = directions(id, d);
	double a = 0;
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		a += direction[d] * direction[d];

	for (unsigned int j = 0; j < numOverlaps; ++j)
	{
		// ray-sphere intersection continues
		double b = intersectionArray->x;
		double c = intersectionArray->y;

		intersectionArray += NUM_BLOCKS * NUM_THREADS;

		// use double precision. 
		// single presion is not enough here.
		double distSqrt = sqrt(b * b - 4 * a * c);
		double q;
		if (b < 0)
			q = (-b - distSqrt) / 2;
		else
			q = (-b + distSqrt) / 2;

		// compute t0 and t1
		DATA_TYPE x0 = (DATA_TYPE)(q / a);
		DATA_TYPE x1 = (DATA_TYPE)(c / q);
		if (x0 > x1)
		{
			// if t0 is bigger than t1 swap them around
			DATA_TYPE temp = x0;
			x0 = x1;
			x1 = temp;
		}

		short distance = interval.numIntervals;
		short low, up, tempIt;
		low = GPUDart::GPUInterval::lower_bound(interval.intervalBegins, interval.next, interval.head, interval.tail, x0, &distance);
		tempIt = low;
		if (low != interval.head)
		{
			tempIt = interval.prev[low];
			distance += 1;
		}
		up = GPUDart::GPUInterval::upper_bound(interval.intervalEnds, interval.next, tempIt, interval.tail, x1, &distance);

		DATA_TYPE min = interval.intervalBegins[interval.head];
		DATA_TYPE max = interval.intervalEnds[interval.prev[interval.tail]];

		tempIt = up;
		// split an interval
		if (up != interval.tail && interval.next[tempIt] == low)
		{
			DATA_TYPE splitEnd = interval.intervalEnds[up];
			interval.intervalEnds[up] = x0;
			interval.insert(low, x1, splitEnd);
		}
		else
		{
			if (x0 > min)
			{
				tempIt = interval.prev[low];
				if (x0 < interval.intervalEnds[tempIt])
					interval.intervalEnds[tempIt] = x0;
			}
			if (x1 < max)
			{
				if (x1 > interval.intervalBegins[up])
					interval.intervalBegins[up] = x1;
			}
			interval.erase(low, up);
		}
	} 
}

__device__ DATA_TYPE raisePower(DATA_TYPE value)
{
	bool invert = (value < 0);
	if (invert) 
		value = -value;
	value = pow(value, (DATA_TYPE)DIMENSIONS);
	if (invert)
		value = -value;
	return value;
}

__device__ DATA_TYPE lowerPower(DATA_TYPE value)
{
	bool invert = (value < 0);
	if (invert) 
		value = -value;
	value = pow(value, c_invNumDimensions);
	if (invert)
		value = -value;
	return value;
}

__global__ void pickPointKernel(unsigned int* pNumFronts, curandState* randStates, GPUSpokeDart::GPUInterval* intervals, 
		GPUDart::GPUMDVector refPoints, GPUDart::GPUMDVector directions, GPUDart::GPUMDVector computedSamples)
{
	unsigned int numFronts = *pNumFronts;
	if (numFronts == 0)
		return;
	unsigned int id = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int frontId = id - (id / numFronts) * numFronts;

	GPUDart::GPUInterval& interval = intervals[id];
	DATA_TYPE refPoint[DIMENSIONS];
	DATA_TYPE direction[DIMENSIONS];

	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		refPoint[d] = refPoints(frontId, d);
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		direction[d] = directions(id, d);

	short numIntervals = interval.numIntervals;
	d_hasSolution[id] = (numIntervals != 0) ? 1 : 0;
	if (numIntervals != 0)
	{
		for (int i = interval.head; i != interval.tail; i = interval.next[i])
		{
			interval.intervalBegins[i] = raisePower(interval.intervalBegins[i]);
			interval.intervalEnds[i] = raisePower(interval.intervalEnds[i]);
		}

		DATA_TYPE totalLength = 0;
		for (int i = interval.head; i != interval.tail; i = interval.next[i])
		{
			totalLength += interval.intervalEnds[i] - interval.intervalBegins[i];
		}

		curandState* localState = randStates + id;
		DATA_TYPE randPoint = RandUtilGpu::randRange(localState, totalLength);

		DATA_TYPE computedSample[DIMENSIONS];
		for (int i = interval.head; i != interval.tail; i = interval.next[i])
		{
			DATA_TYPE length = interval.intervalEnds[i] - interval.intervalBegins[i];
			if (randPoint <= length)
			{
				DATA_TYPE t = lowerPower(interval.intervalBegins[i] + randPoint);
				for (unsigned int d = 0; d < DIMENSIONS; ++d)
					computedSample[d] = refPoint[d] + t * direction[d];

				//printf("[%d] t : %f (%f %f) = (%f %f) * %f (%f %f)\n", id, t, computedSample[0], computedSample[1], refPoint[0], refPoint[1], t, direction[0], direction[1]);
				break;
			}
			randPoint -= length;
		}

		for (unsigned int d = 0; d < DIMENSIONS; ++d)
			computedSamples(id, d) = computedSample[d];

		//printf("[%d] picked (%f, %f)\n", id, computedSample[0], computedSample[1]);
	}
}

__global__ void addSamplesKernel(unsigned int* pNumFronts, GPUDart::GPUMDVector samples, 
		unsigned int *pNumSamples, unsigned int* pFailureCounts, GPUDart::GPUMDVector computedSamples,
		unsigned int* d_gridcells, unsigned int* d_numGridCellElements)
{
	if (*pNumFronts == 0)
		return;
	unsigned int id = threadIdx.x + blockIdx.x * blockDim.x;

	__shared__ unsigned char sharedHasSample[NUM_THREADS];
	__shared__ VectorGPU sharedSamples[NUM_THREADS];

	sharedHasSample[threadIdx.x] = d_hasSolution[id];
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
		sharedSamples[threadIdx.x][d] = computedSamples(id, d);

	__syncthreads();

	VectorGPU diff;

	// discard in block
	for (int i = 0; i < blockDim.x; ++i)
	{
		if (i < threadIdx.x && sharedHasSample[i] && sharedHasSample[threadIdx.x])
		{
			VectorGPU otherSample = sharedSamples[i];
			VectorGPU thisSample = sharedSamples[threadIdx.x];
			for (int d = 0; d < DIMENSIONS; ++d)
				diff[d] = otherSample[d] - thisSample[d];
			if (diff.getSquaredNorm() < c_sqRadius)
			{
				sharedHasSample[threadIdx.x] = false;
				//printf("[%d] discarded\n", threadIdx.x);
				break;
			}
		}
	}

	__syncthreads();

	if (threadIdx.x == 0)
	{
		unsigned int numCheckedSamples = d_checkedSamples[blockIdx.x];
		unsigned int numBlockSamples = 0;
		unsigned int blockSampleIndices[NUM_THREADS];
		for (unsigned int s = 0; s < NUM_THREADS; ++s)
		{
			if (sharedHasSample[s])
			{
				blockSampleIndices[numBlockSamples++] = s;
			}
		}

		unsigned int numSamples;
		unsigned int numUncheckedSamples;

		// acquire lock
		while (atomicExch(&d_lock, 1) != 0)
		{
		}
		numSamples = *pNumSamples;
		numUncheckedSamples = numSamples;

		if (numUncheckedSamples != numCheckedSamples)
		{
			for (unsigned int i = numCheckedSamples; i < numUncheckedSamples; ++i)
			{
				VectorGPU sample;
				for (unsigned int d = 0; d < DIMENSIONS; ++d)
					sample[d] = samples(i, d);
				for (unsigned int s = 0; s < numBlockSamples; )
				{
					unsigned int sampleIndex = blockSampleIndices[s];
					VectorGPU refPoint = sharedSamples[sampleIndex];

					for (unsigned int d = 0; d < DIMENSIONS; ++d)
						diff[d] = sample[d] - refPoint[d];

					if (diff.getSquaredNorm() < c_sqRadius)
					{
						if (s != numBlockSamples - 1)
							blockSampleIndices[s] = blockSampleIndices[numBlockSamples - 1];
						--numBlockSamples;
					}
					else
					{
						++s;
					}
				}
			}
			numCheckedSamples = numUncheckedSamples;
		}

		for (unsigned int s = 0; s < numBlockSamples; ++s)
		{
			unsigned int sampleIndex = blockSampleIndices[s];
			VectorGPU newSample = sharedSamples[sampleIndex];
			for (unsigned int d = 0; d < DIMENSIONS; ++d)
				samples(numSamples, d) = newSample[d];
#ifdef USE_GRID_BASED_NEIGHBOR_SEARCH
				addToGrid(newSample, numSamples, d_gridcells, d_numGridCellElements);
#endif
				//printf("[%d] adds (%f, %f)as sample %d\n", sampleIndex, newSample[0], newSample[1], numSamples);
			++numSamples;
		}

		(*pNumSamples) = numSamples;

		// release lock
		atomicExch(&d_lock, 0);
	}
}


__global__ void updateFront(unsigned int* pNumFronts, unsigned int *pNumSamples, unsigned int* pFailureCounts, unsigned int* pNextFront, unsigned int numBlocks)
{
	unsigned int frontId = threadIdx.x;

	unsigned int numFronts = *pNumFronts;
	unsigned int numSamples = *pNumSamples;

	if (numFronts == 0)
		return;

	__shared__ unsigned int sharedFronts[NUM_BLOCKS * NUM_THREADS];

	// update failure count & advance front
	for (; frontId < numFronts; frontId += blockDim.x)
	{
		sharedFronts[frontId] = d_front[frontId];
		unsigned int failureCount = pFailureCounts[frontId];
		// unroll?
		for (unsigned int i = frontId; i < numBlocks * NUM_THREADS; i += numFronts)
		{
			unsigned char trialHasSolution = d_hasSolution[i];
			if (trialHasSolution)
				failureCount = 0;
			else
				++failureCount;
		}

		if (failureCount >= c_m)
		{
			sharedFronts[frontId] = atomicInc(pNextFront, INT_MAX);
			failureCount = 0;

			//printf("replaced front [%d] : %d->%d\n", frontId, old, sharedFronts[frontId]);
		}
		pFailureCounts[frontId] = failureCount;
	}
	__syncthreads();

	// remove empty fronts
	unsigned int backup;
	if (threadIdx.x == 0)
	{
		for (unsigned int i = 0; i < numFronts; ++i)
		{
			if (sharedFronts[i] >= numSamples)
			{
				--numFronts;
				if (i != numFronts)
				{
					//printf("front [%d] removed : %d\n", i, sharedFronts[i]);
					sharedFronts[i] = sharedFronts[numFronts];
					pFailureCounts[i] = pFailureCounts[numFronts];
					--i;
				}
				pFailureCounts[numFronts] = 0;
			}
		}
		backup = sharedFronts[0];
		sharedFronts[0] = numFronts;
	}
	__syncthreads();
	if (threadIdx.x != 0)
		numFronts = sharedFronts[0];
	__syncthreads();
	if (threadIdx.x == 0)
		sharedFronts[0] = backup;

	// add new fronts
	unsigned int maxFront = numBlocks * NUM_THREADS;
	unsigned int nextFront = min(*pNextFront, numSamples);
	unsigned int numAddedFronts = min(numSamples - nextFront, maxFront - numFronts);
	for (frontId = numFronts + threadIdx.x; frontId < numFronts + numAddedFronts; frontId += blockDim.x)
	{
		sharedFronts[frontId] = nextFront + frontId - numFronts;
	}
	__syncthreads();

	numFronts += numAddedFronts;

	// write
	for (frontId = threadIdx.x; frontId < numFronts; frontId += blockDim.x)
	{
		d_front[frontId] = sharedFronts[frontId];
	}

	if (threadIdx.x == 0)
	{
		*pNumFronts = numFronts;
		*pNextFront = nextFront + numAddedFronts;
	}
}

};

////////////////////////////////////////////////////////////////////////////////

GPUSpokeDart::GPUSpokeDart()
{

}

GPUSpokeDart::~GPUSpokeDart()
{
	CPUFREE(h_buffer_);

	d_samples_.destroy();
	d_refPoints_.destroy();
	d_directions_.destroy();
	d_computedSamples_.destroy();

	GPUFREE(d_randStates_);
	GPUFREE(d_intervals_);
	GPUFREE(d_numSamples_);
	GPUFREE(d_failureCount_);

#ifdef USE_GRID_BASED_NEIGHBOR_SEARCH
	GPUFREE(GPUSpokeDartCuda::d_numGridCellElements);
	GPUFREE(GPUSpokeDartCuda::d_gridcells);
#endif
}

bool GPUSpokeDart::init(unsigned int m, DATA_TYPE radius, unsigned int d)
{
	if (d != DIMENSIONS)
	{
		printf("Dimension mismatch!!!\n");
		return false;
	}
#ifdef USE_GRID_BASED_NEIGHBOR_SEARCH
	if (d != 2)
	{
		printf("Grid-based neighbor search only can be used in 2D\n");
		return false;
	}
#endif
		

	cudaDeviceSetSharedMemConfig(cudaSharedMemBankSizeEightByte);

	h_m_ = m;
	h_radius_ = radius;
	h_dimensions_ = d;

	d_samples_.init(MAX_SAMPLES, h_dimensions_);
	d_refPoints_.init(MAX_FRONT, h_dimensions_);
	d_directions_.init(NUM_THREADS * NUM_BLOCKS, h_dimensions_);
	d_computedSamples_.init(NUM_THREADS * NUM_BLOCKS, h_dimensions_);

	TOGPU_CONSTANT(GPUSpokeDartCuda::c_m, &m, sizeof(unsigned int));
	TOGPU_CONSTANT(GPUSpokeDartCuda::c_radius, &radius, sizeof(DATA_TYPE));
	DATA_TYPE sqRadius = radius * radius;
	TOGPU_CONSTANT(GPUSpokeDartCuda::c_sqRadius, &sqRadius, sizeof(DATA_TYPE));
	DATA_TYPE invDimensions = (DATA_TYPE)1 / DIMENSIONS;
	TOGPU_CONSTANT(GPUSpokeDartCuda::c_invNumDimensions, &invDimensions, sizeof(DATA_TYPE));

	GPUMALLOC(&d_randStates_, sizeof(curandState) * NUM_THREADS * NUM_BLOCKS);
	GPUMALLOC(&d_intervals_, sizeof(GPUInterval) * NUM_THREADS * NUM_BLOCKS);
	GPUMALLOC(&d_numSamples_, sizeof(unsigned int) * 3);
	d_numFront_ = (unsigned int*)d_numSamples_ + 1;
	d_nextFront_ = (unsigned int*)d_numSamples_ + 2;
	GPUMALLOC(&d_failureCount_, sizeof(unsigned int) * MAX_FRONT);

	CPUMALLOC(&h_buffer_, sizeof(unsigned int) * 3);

	// add the first sample
	VectorGPU firstSample;
	VectorGPU* d_firstSample;
	GPUMALLOC(&d_firstSample, sizeof(VectorGPU));
	for (unsigned int d = 0; d < DIMENSIONS; ++d)
	{
		firstSample[d] = RandUtil::rand0To1();
	}
	TOGPU(d_firstSample, &firstSample, sizeof(VectorGPU));

	h_numSamples_ = 1;
	TOGPU((void*)d_numSamples_, &h_numSamples_, sizeof(unsigned int));

	h_numFront_ = 1;
	TOGPU((void*)d_numFront_, &h_numFront_, sizeof(unsigned int));
	GPUMEMSET(d_failureCount_, 0, sizeof(unsigned int) * MAX_FRONT);

#ifdef USE_GRID_BASED_NEIGHBOR_SEARCH
	{
		// set grid (2D)
		DATA_TYPE h_cellSize = 1.5 * h_radius_;
		int h_gridMax = (int)ceil((DATA_TYPE)1 / h_cellSize);
		TOGPU_CONSTANT(GPUSpokeDartCuda::d_cellSize, &h_cellSize, sizeof(DATA_TYPE));
		TOGPU_CONSTANT(GPUSpokeDartCuda::d_gridMax, &h_gridMax, sizeof(int));
		GPUMALLOC(&GPUSpokeDartCuda::d_numGridCellElements, sizeof(unsigned int) * h_gridMax * h_gridMax);
		GPUMEMSET(GPUSpokeDartCuda::d_numGridCellElements, 0, sizeof(unsigned int) * h_gridMax * h_gridMax);
		GPUMALLOC(&GPUSpokeDartCuda::d_gridcells, sizeof(unsigned int) * h_gridMax * h_gridMax * MAX_GRID_CELL_ELEMENTS);

		// set first value
		int x = (int)(firstSample[0] / h_cellSize);
		int y = (int)(firstSample[1] / h_cellSize);
		unsigned int count = 1;
		unsigned int sampleIndex = 0;
		TOGPU(GPUSpokeDartCuda::d_numGridCellElements + (y * h_gridMax + x), &count, sizeof(unsigned int));
		TOGPU(GPUSpokeDartCuda::d_gridcells + (y * h_gridMax + x) * MAX_GRID_CELL_ELEMENTS, &sampleIndex, sizeof(unsigned int));
	}
#endif


	GPUSpokeDartCuda::initRandomStates <<< NUM_BLOCKS, NUM_THREADS >>> (d_randStates_, d_nextFront_, d_samples_, d_firstSample);

	CUDA_CHECK_ERROR();

	GPUFREE(d_firstSample);

	return true;
}

void GPUSpokeDart::setupDarts()
{
	int numBlocks = useMultipleBlocks_ ? NUM_BLOCKS : 1;
	GPUSpokeDartCuda::setupGPUDarts <<< numBlocks, NUM_THREADS >>> (d_numFront_, d_samples_, d_randStates_, d_refPoints_, d_directions_, d_intervals_);

	CUDA_CHECK_ERROR();
}

void GPUSpokeDart::updateIntervals()
{
	int numBlocks = useMultipleBlocks_ ? NUM_BLOCKS : 1;
	GPUSpokeDartCuda::updateIntervalsKernel <<< numBlocks, NUM_THREADS >>> (d_numFront_, d_samples_, d_intervals_, d_numSamples_,
		d_refPoints_, d_directions_,
		GPUSpokeDartCuda::d_gridcells, GPUSpokeDartCuda::d_numGridCellElements);

	CUDA_CHECK_ERROR();
}

void GPUSpokeDart::updateIntervals2()
{
	int numBlocks = useMultipleBlocks_ ? NUM_BLOCKS : 1;
	GPUSpokeDartCuda::updateIntervals2Kernel <<< numBlocks, NUM_THREADS >>> (d_numFront_, d_samples_, d_intervals_, d_numSamples_,
		d_refPoints_, d_directions_,
		GPUSpokeDartCuda::d_gridcells, GPUSpokeDartCuda::d_numGridCellElements);

	CUDA_CHECK_ERROR();
}

void GPUSpokeDart::pickPoint()
{
	int numBlocks = useMultipleBlocks_ ? NUM_BLOCKS : 1;
	GPUSpokeDartCuda::pickPointKernel <<< numBlocks, NUM_THREADS >>> (d_numFront_, d_randStates_, d_intervals_, 
		d_refPoints_, d_directions_, d_computedSamples_);

	CUDA_CHECK_ERROR();
}

void GPUSpokeDart::addSamples()
{
	static int doSync = 100;

	if (useMultipleBlocks_ && doSync > 0)
		cudaDeviceSynchronize();

	int numBlocks = useMultipleBlocks_ ? NUM_BLOCKS : 1;
	GPUSpokeDartCuda::addSamplesKernel <<< numBlocks, NUM_THREADS >>> (d_numFront_, d_samples_, d_numSamples_,
		d_failureCount_, d_computedSamples_,
		GPUSpokeDartCuda::d_gridcells, GPUSpokeDartCuda::d_numGridCellElements);

	if (useMultipleBlocks_ && doSync > 0)
	{
		cudaDeviceSynchronize();
		--doSync;
	}

	CUDA_CHECK_ERROR();

	GPUSpokeDartCuda::updateFront <<< 1, NUM_THREADS >>> (d_numFront_, d_numSamples_, d_failureCount_, d_nextFront_, numBlocks);

	CUDA_CHECK_ERROR();
}

bool GPUSpokeDart::terminationCheck()
{
	// TODO:
	return h_numFront_ == 0;
}

unsigned int GPUSpokeDart::getNumGeneratedSamples()
{
	FROMGPU(h_buffer_, d_numSamples_, sizeof(unsigned int) * 3);
	h_numSamples_ = h_buffer_[0];
	h_numFront_ = h_buffer_[1];
	h_numNextFront_ = h_buffer_[2];
	//printf("# samples : %d # front : %d\n", h_numSamples_, h_numFront_);

	return h_numSamples_;
}

void GPUSpokeDart::downloadSamples()
{
	FROMGPU(h_buffer_, (unsigned int*)d_numSamples_, sizeof(unsigned int));
	h_numSamples_ = h_buffer_[0];

	if (h_numSamples_ == 0)
		return;

	DATA_TYPE* h_temp = new DATA_TYPE[d_samples_.size_ * h_dimensions_];
	FROMGPU(h_temp, d_samples_.p_, sizeof(DATA_TYPE) * d_samples_.size_ * h_dimensions_);
	VectorGPU* h_samples = new VectorGPU[h_numSamples_];
	for (unsigned int i = 0; i < h_numSamples_; ++i)
	{
		for (unsigned int d = 0; d < h_dimensions_; ++d)
		{
			h_samples[i][d] = h_temp[d * d_samples_.size_ + i];
		}
	}

	/*
	for (int i = 0; i < h_numSamples_; ++i)
	{
		printf("[%d] ", i);
		for (int d = 0; d < h_dimensions_; ++d)
		{
			printf("%f ", h_samples[i].elements_[d]);
		}
		printf("\n");
	}
	*/

	if (validateResults_)
	{
		printf("validation\n");
		for (unsigned int i = 0; i < h_numSamples_; ++i)
		{
			for (unsigned int j = i + 1; j < h_numSamples_; ++j)
			{
				DATA_TYPE sq = (h_samples[i] - h_samples[j]).getSquaredNorm();
				if (sq < h_radius_ * h_radius_ - DATA_TYPE_EPS)
				{
					printf("!!!invalid %d-%d dist : %f\n", i, j, sqrt(sq));
				}
			}
		}
	}

	printf("# of samples : %d\n", h_numSamples_);

	delete[] h_samples;
	delete[] h_temp;
}

////////////////////////////////////////////////////////////////////////////////

